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ABSTRACT 

Following the optical imaging of exoplanet candidate Fomalhaut b (Fom b) , we present a numerical 
model of how Fomalhaut’s debris disk is gravitationally shaped by a single interior planet. The model 
is simple, adaptable to other debris disks, and can be extended to accommodate multiple planets. 
If Fom b is the dominant perturber of the belt, then to produce the observed disk morphology it 
must have a mass M p \ < 3Mj, an orbital semimajor axis a p i > 101.5 AU, and an orbital eccentricity 
e p i = 0.11-0.13. These conclusions are independent of Fom b’s photometry. To not disrupt the 
disk, a greater mass for Fom b demands a smaller orbit farther removed from the disk; thus, future 
astrometric measurement of Fom b’s orbit, combined with our model of planet-disk interaction, can 
be used to determine the mass more precisely. The inner edge of the debris disk at a ~ 133 AU lies 
at the periphery of Fom b’s chaotic zone, and the mean disk eccentricity of e ss 0.11 is secularly 
forced by the planet, supporting predictions made prior to the discovery of Fom b. However, previous 
mass constraints based on disk morphology rely on several oversimplifications. We explain why our 
constraint is more reliable. It is based on a global model of the disk that is not restricted to the 
planet’s chaotic zone boundary. Moreover, we screen disk parent bodies for dynamical stability over 
the system age of ~100 Myr, and model them separately from their dust grain progeny; the latter’s 
orbits are strongly affected by radiation pressure and their lifetimes are limited to ~0.1 Myr by 
destructive grain-grain collisions. The single planet model predicts that planet and disk orbits be 
apsidally aligned. Fomalhaut b’s nominal space velocity does not bear this out, but the astrometric 
uncertainties are difficult to quantify. Even if the apsidal misalignment proves real, our calculated 
upper mass limit of 3 Mj still holds. Parent bodies are evacuated from mean-motion resonances with 
Fom b; these empty resonances are akin to the Kirkwood gaps opened by Jupiter. The belt contains 
at least 3 M® of solids that are grinding down to dust, their velocity dispersions stirred so strongly 
by Fom b that collisions are destructive. Such a large mass in solids is consistent with Fom b having 
formed in situ. 

Subject headings: stars: planetary systems — stars: circumstellar matter - - planetary systems: pro- 
toplanetary disks — celestial mechanics — stars: individual (Fomalhaut) 


1. INTRODUCTION 

A common proper motion companion to Fomalhaut 
has been imaged by Kalas et al. (2008, hereafter K08) 
using the Hubble Space Telescope Advanced Camera for 
Surveys (HST ACS) coronagraph. Fomalhaut b (Fom 
b) orbits interior to the system’s well-known circumstel- 
lar belt of dust (e.g., Holland et al. 1998; Kalas et al. 
2005, hereafter K05, and references therein). While Fom 
b’s ultra-low luminosity leaves little doubt that it is of 
remarkably low mass, sitting well below the regime of 
brown dwarfs, the question remains: Just how low is its 
mass? 

Based on the observed broadband spectrum of Fom b, 
K08 estimate an upper limit on the mass M p \ of about 
3Mj. We recapitulate their reasoning as follows. Fom 
b is detected in HST’s F814W (0.7-0. 9 /irn) and F606W 
(0.45-0.7 jim) passbands in 2006. The F606W flux is 
variable; the flux in 2006 was about half that in 2004. 
Observations in 2005 with Keck in H band (1.5-1. 8 jim) 
and in 2008 with Gemini in L-prime (3.2-4 jim) give only 
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upper limits. 

The F814W flux (observed only in 2006; imaging was 
not attempted in 2004 for this passband) can be repro- 
duced by thermal emission from a 2-4Mj, 200-Myr-olcl 
planet that formed by core accretion and is of supersolar 
metallicity (Hubickyj et al. 2005; Fortney et al. 2008). 
Unfortunately, this same model atmosphere underpre- 
dicts, by more than an order of magnitude, the F606W 
fluxes. At the same time it overpredicts the H-band 
3cr upper limit by a factor of a few, and is marginally 
consistent with the L-prime 3cr upper limit. From con- 
siderations outlined by K08, we can construct two, not 
entirely exclusive hypotheses. Hypothesis one: F814W 
still traces planetary thermal emission, F606W is con- 
taminated by variable Ha emission, and the atmospheric 
model requires revision in H band. The Ha hypothesis 
is inspired by variable Ha emission from chromospheri- 
cally active and/or accreting low-mass stars and brown 
dwarfs. The puzzling brown dwarf GQ Lup B offers a 
possible precedent (Marois et al. 2007); in either the case 
of Fom b or that of GQ Lup B, the Ha luminosity nec- 
essary to explain the anomalously large F606W flux is 
~1% that of the bolometric luminosity (unfortunately in 
neither case has an optical spectrum been taken) . As for 
uncertainties in model exoplanet atmospheres (Fortney 
et al. 2008; Burrows et al. 2003), these appear greater 
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in H band than in F814W; the various models disagree 
with each other at near-infrared wavelengths by factors 
of a few. Hypothesis two: both F814W and F606W are 
contaminated by starlight reflected off an optically thick 
circumplanetary disk, and the F606W variability arises 
from variable disk accretion onto the planet, possibly im- 
plicating Hcc again. To explain the detected fluxes in 
2006 using reflected light alone, such a disk would have to 
be comparable in size to the orbits of Jupiter’s Galilean 
satellites. It is not clear, however, how such a disk would 
survive for the system age of ~100 Myr; protosatellite 
disks are thought to evolve on timescales shorter than a 
few Myr (e.g., Canup & Ward 2002; Mosqueira & Estrada 
2003). 

Since both hypotheses admit the possibility of addi- 
tional sources of luminosity apart from planetary ther- 
mal emission, the 3Mj inferred from the observed F814W 
flux should be considered an upper limit. As a whole, 
this interpretation of Fom b’s photometry is preliminary, 
subject to significant revision as both observations and 
theory improve. 

An alternative route to probing the properties of Fom 
b is to exploit the morphology of the circumstellar belt. 
Prior to the discovery of Fom b, this approach was taken 
by Quillen (2006, hereafter Q06), who built on earlier 
work by Wyatt et al. (1999). The striking intrinsic el- 
lipticity of the belt of e « 0.11 can be forced by secular 
gravitational interaction with a single planet on a simi- 
larly eccentric orbit interior to the belt. 4 Fomalhaut b 
was predicted by Q06 to occupy an orbit of eccentricity 
0.1 and semimajor axis ~119 AU, and to have a mass 
between about 2 x 10 -5 and 7 x 10 -5 that of the cen- 
tral star. For a stellar mass of M* = 2.3 M 0 , this range 
corresponds to ~0.05-0.2 Mj. 

The predictions of Q06 rest on the idea that the belt in- 
ner edge, at senrinrajor axis dinner = 133 AU, is located at 
the outer boundary of Fom b’s chaotic zone. The chaotic 
zone is a swath of space enclosing the planet’s orbit in 
which test particle orbits are chaotic and short-lived, as a 
consequence of overlapping first-order mean-motion reso- 
nances (Wisdom 1980). The semimajor axes a c h a otic that 
bound this chaotic zone are displaced to either side of the 
planet’s semimajor axis a p i by 

Adchaotic = |dchaotic U p l| ~ 1-3 /X ^ (1) 

where /x = M p i/M*. The coefficient of 1.3 arises from 
Wisdom’s approximate scaling theory. Though (1) is de- 
rived for the case of a planet that occupies a circular 
orbit and that interacts with test particles on nearly cir- 
cular orbits, Quillen & Faber (2006) find that the /x 2 / 7 
scaling law holds also for a planet on a moderately eccen- 
tric orbit, interacting with particles on secularly forced 
eccentric orbits. Quillen & Faber (2006) prefer, how- 
ever, the coefficient of 1.5 that originates from the nu- 
merical integrations and eccentricity growth criterion of 
Duncan et al. (1989). 5 * Our work supports the assign- 
ment dinner = d c h a otic made by Q06. However, we will 

4 See the textbook by Murray & Dermott (2000) for an introduc- 
tion to secular perturbation theory, which essentially treats masses 
as orbit-averaged elliptical wires. 

5 Duncan et al. (1989) also provide a simple explanation of be- 

havior within the chaotic zone (for the case of circular orbits): 
particles residing there are perturbed so strongly by the planet 

at each conjunction that successive conjunctions occur at uncor- 


calculate still a third coefficient that best reproduces the 
HST images of Fomalhaut’s belt; see our equation (16). 

The detection of Fom b orbiting just interior to the 
dust belt appears to confirm the expectation of K05, 
made quantitative by Q06, that a planet is responsible 
for truncating the inner edge of the belt and endowing 
the belt with its eccentric shape. However, a closer ex- 
amination reveals a potential problem with the hypothe- 
sis that Fom b is solely responsible for the observed belt 
morphology. This picture predicts that the planet’s orbit 
be apsidally aligned and coplanar with (perfectly nested 
inside) that of the belt. If that were the case, Fom b’s 
observed position would place it somewhat past orbital 
quadrature, near a true anomaly of 109 deg. The prob- 
lem is that the expected orbital velocity near quadrature 
(about 4.2 krn/s for a stellar mass of 2.3 Mq and a semi- 
major axis of 115 AU) is lower than K08’s estimate for 
the actual, deprojected space velocity of Fom b (5.5t 0 'y 
km/s). The large observed space velocity suggests that 
Fom b is currently located closer to periastron, and that 
planet and belt orbits are not apsidally aligned. 

Should the apsidal misalignment be taken seriously and 
the single planet picture abandoned? It is hard to say. 
With Fom b detected at only two epochs so far, it is im- 
possible to derive a unique orbit. The systematic errors 
underlying the measured orbital velocity — in particular 
frame registration errors that arise because the star is 
hidden behind the coronagraphic spot — are difficult to 
quantify with confidence. More precise statements about 
Fom b’s orbit will have to await future astrometry. Our 
current prejudice is to say that the observational uncer- 
tainties are large, and that apsidal alignment remains a 
possibility. Given the observed proximity of Fom b to the 
belt, it would be strange if Fom b were not the dominant 
perturber (see also the last paragraph of §4.2). 

The unresolved issue of apsidal alignment notwith- 
standing, we can still place an upper limit to the mass 
of Fom b by assuming that it is solely responsible for 
shaping the belt. Broadly speaking, more planets render 
more of orbital phase space chaotic. Therefore Fom b, if 
abetted by other planets, could have a smaller mass than 
that derived under the single planet assumption and still 
truncate the belt at its observed inner edge. The mass 
derived under the single planet assumption would be an 
upper limit. 

Within the single planet picture, the published dynam- 
ical constraint on the planet’s mass is 0.05 < M p i(Mj) < 
0.2 (Q06). But there are several reasons to question this 
result: 

1. Belt properties depend weakly on planet mass, and 
therefore uncertainties in the former are magni- 
fied when constraining the latter. For example, 
Aa c haotici which governs the relative locations of 
the belt inner edge and the planet, scales only as 
MJ . Similarly weak powers describe how the ve- 
locity dispersion at the boundary of the chaotic 
zone depends on M p \; Q06 parlays this dependence 
into an upper limit on M p i, on the grounds that 
too large a velocity dispersion violates the observed 

related longitudes; consequently the particle undergoes a random 
walk in semimajor axis and eccentricity, with steps in the walk 
corresponding to impulses at every conjunction. 
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sharpness of the belt’s inner edge. We will see in 
our work that this argument does not fully capture 
the actual behavior because it neglects how belt 
particles located some distance from the chaotic 
zone boundary influence the observed sharpness of 
the edge. In other words, sharpness can only be 
reliably computed with a global model, not one 
that examines the chaotic zone boundary exclu- 
sively. See Figure 8 and the related discussion in 
§3.3.1. 

2. The published dynamical upper mass limit is based 
on the purely gravitational, collisionless behavior of 
test particles. But the HST observations of Fomal- 
haut’s belt are of dust grains, which are influenced 
by stellar radiation pressure and interparticle col- 
lisions. As explained in §2 (see also Strubbe & 
Cliiang 2006), the scattered light observations are 
likely dominated by the smallest grains still bound 
to the star. Of all the solid material orbiting Foma- 
lhaut, such grains have the largest ratios of surface 
area to mass and are the most seriously affected by 
radiative forces and collisions. 

3. The lower dynamical mass limit is also suspect be- 
cause it is premised on particles colliding indestruc- 
tibly and diffusing as members of a viscous circular 
ring. The argument underlying the lower limit is 
that the planet mass cannot be too small lest dust 
grains diffuse into the chaotic zone by interparti- 
cle collisions, before the planet can gravitationally 
eject them (see also Quillen 2007). But dust par- 
ticles, colliding at minimum speeds of ~100 m/s 
(see our §2), likely shatter one another. Moreover, 
their trajectories are highly elliptical as a conse- 
quence of radiation pressure (Strubbe & Chiang 
2006). Their dynamics seem poorly described by 
a diffusion equation that conserves particle num- 
ber, has constant diffusivity, and assumes circular 
orbits. 

In this work we present a new numerical model of the 
Fomalhaut dust belt under the single planet assumption. 
It accounts not only for gravitational sculpting by Fom 
b, but also for the finite collisional lifetime of dust, the 
size distribution of grains, and radiation forces (includ- 
ing Poynting-Robertson drag, though it is of minor im- 
portance compared to other perturbations). Though we 
do not go as far as generating a model scattered light 
image to compare with observations, we take the first 
step towards this goal by calculating the detailed shape 
of the belt’s vertical optical depth, rj_, as a function of 
semimajor axis a. We compare this optical depth pro- 
file with the corresponding profile of the K05 scattered 
light model, which in turn was fitted directly to the HST 
images. From this comparison we constrain the mass 
and orbit of Fom b: we calculate the possible combina- 
tions of planet mass M p \, orbital semimajor axis a p i, and 
orbital eccentricity e p i. These results are independent of 
any measurement of Fom b, in particular of its spectrum. 
Only as a final extra step do we use information on Fom 
b’s observed stellocentric distance to see if some of our 
mass-orbit combinations might be ruled out. 

Our work corroborates the overall picture of Q06 — 
that the planet’s chaotic zone truncates the inner edge 


of the belt, and that the planet’s orbital eccentricity in- 
duces a similar mean eccentricity in the belt (Wyatt et al. 
1999)- but we introduce enough improvements, espe- 
cially with regards to our separate handling of unobserv- 
able parent bodies and observable dust grains, that our 
dynamical constraint on Fom b’s mass supersedes that 
of Q06. Moreover, if future astrometry reveals that that 
there is no significant apsidal misalignment, then the pre- 
cise relationship we derive between Fom b’s mass and its 
orbit can be used to determine the former using the lat- 
ter. Our model is simple, robust to uncertainties in input 
parameters, and applicable to other systems. Though 
we restrict consideration to a single interior planet, ad- 
ditional perturbers can be added. 

In §2 we present an overview of the Fomalhaut belt, 
estimating its physical properties to order-of-magnitude 
accuracy. These estimates inform our choices for the in- 
put parameters of our numerical model; that model, and 
its output, are detailed in §3. We summarize and discuss 
our results — describing also the curious anomalous accel- 
eration of Fomalhaut measured by the Hipparcos satel- 
lite, and how other planetary companions in addition to 
Fom b affect our conclusions — in §4. 

2. ORDERS OF MAGNITUDE 

We derive basic properties of the Fomalhaut dust belt, 
working as much as possible from first principles and 
direct observations. The conclusions reached in the fol- 
lowing subsections form the basis of our numerical model 
in §3. 

2.1. Bound dust grains present absorbing, geometric 
cross sections 

Fomalhaut is a spectral type A star of luminosity 
L* = 16 Lq, mass M* « 2.3 M 0 , and age t age = 200± 100 
Myr (Barrado y Navascues et al. 1997; Barrado y Navas- 
cues 1998). The circumstellar debris emits a fractional 
infrared excess of Tir/T* = 4.6 x 10 -5 (Song et al. 2001). 
K05 report that from 0.6 to 0.8 ^m wavelengths, the 
ring has a spatially integrated apparent magnitude of 
7U a pp = 16.2, nearly entirely due to reflected starlight; 
for the star, m app ,* = 1.12; therefore the fractional re- 
flected luminosity is L re f/T* « 10” 6 . If we ignore order 
unity effects introduced by anisotropy in the scattering 
phase function (see K05 for a fit to this phase function), 
the dust albedo is of order L re f/LiR « 0.02 — the grains 
are nearly purely absorbing. 

Dust is generated by the collisional comminution of 
larger parent bodies. The largest parent bodies sit at the 
top of the collisional cascade and have lifetimes against 
collisional disruption equal to t age . We take these largest 
parents to occupy a torus of radius R « 140 AU, annular 
width A R < R, and uniform vertical thickness 2 H < R, 
to match the observations of K05. Strubbe & Chiang 
(2006) refer to this torus as the “birth ring.” In principle, 
a grain of given size that is born into the torus is lost from 
it in one of four ways: collisional comminution, expulsion 
by radiation pressure, ejection by planetary scattering, 
or orbital decay by Poynting-Robertson (PR) drag. In 
practice, for many debris disks, the first two channels are 
more important than the fourth (Wyatt 2005; Strubbe 
& Chiang 2006; Thebault & Wu 2008). The ratio of the 
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force of radiation pressure to that of stellar gravity is 

f- 3L 


I67 vcGM^ps 


( 2 ) 


for a geometrically absorbing grain of internal density 
p ss 1 gcnr -3 and radius s, where c is the speed of light 
and G is the gravitational constant. Grains are unbound 
from the star when /3 > 1/2, 6 i.e. , when 


S ^ ^blow 


3 L* 


87 rcGM*p 


= 8 fim. . 


(3) 


Because Sbiow is much greater than the submicron wave- 
lengths at which the star principally emits, cross sections 
for absorption of radiation by bound grains are indeed 
practically geometric, as (2) assumes. 


2.2. Radiation Pressure Delivers Grains onto Eccentric, 
Long-Period Orbits 

Upon release from parent bodies moving on circular 
orbits, bound grains move on orbits of eccentricity 

ft Sbiow / 

1 — /3 2s — Sbiow 

and semimajor axis 


These expressions, which serve only to guide and which 
are not used in our more precise numerical models, ignore 
the parent-grain relative velocities with which grains are 
born, and also any eccentricity in the parent body or- 
bit. Neither of these errors is serious for the large grain 
eccentricities of interest here (e > 1/3; see §2.5). The 
main point is that radiation pressure flings smaller bound 
grains born in the torus onto more eccentric, longer pe- 
riod orbits. 


2.3. Collisions Between Grains Are Destructive 

Colliding belt particles will chip and shatter one an- 
other. For an angular orbital velocity LIr at semi-major 
axis R , the relative velocity between grains is at least as 
large as the vertical velocity dispersion, ~ HQr = 100 
m/s. To place this velocity into some perspective, we 
note that it is comparable to the maximum flow speeds 
of commercial sandblasting machines. 7 

As a consequence of radiation pressure, many of the 
grains will be travelling on bound orbits having eccen- 
tricities on the order of unity (§2.2). Accounting for 
orbital eccentricities (in-plane velocity dispersion) only 
increases collision velocities, up to a maximum given by 
the local Kepler speed of RDr = 4 krn/s. This maxi- 
mum is comparable to elastic wave speeds in rock and 
will result in catastrophic shattering. 

6 The critical value of /3 = 1/2 applies strictly to dust grains re- 
leased from parent bodies that move initially on circular orbits. For 
dust grains released from parent bodies moving on mildly eccentric 
orbits, as is the case in Fomalhaut’s eccentric belt, the critical 3 
varies with the longitude of release, but is still near 1/2. 

7 See, e.g., http://www.nauticexpo.com/boat- 

manufacturer / sandblasting-machine- 199 1 1 . html. 


2.4. Optical Depths: Radial and Vertical 

Since the grains present largely geometric cross sec- 
tions for absorption of starlight, Lir./L* equals the frac- 
tion of the celestial sphere, centered on the star, that is 
subtended by dust grains: 

L m 2tt R x 2 H xt r H 

TT = irf =r tb ’ (6) 

where tr <C 1 is the radial geometric optical depth 
through the torus. K05 give a model-dependent aspect 
ratio of H/R = 0.025; then tr = 1.8 x 10 -3 . 

The vertical optical depth (measured perpendicular to 
the belt midplane) is 


2 H L m 2 R 

T± ~ Tr a r ~ “a r 


(7) 


independent of H. Again from the scattered light obser- 
vations, A R/R « 0.17 (K05), whence r± = 5.4 x 10“ 4 . 


2.5. Observable Grains in Fomalhaut’s Belt are Bound, 
and Their Lifetime is Set by Collisions, Not by PR 
Drag 

Unbound grains, for which /? > 1/2, exit the torus after 
a local dynamical time, 

tesc($ ^ Sbiow) ^ 77 ^ 200 yr . (8) 

' ‘R 

Because the lifetime of unbound grains, t esc , is much 
shorter than the lifetime of bound grains — the latter life- 
time is set by collisional disruption; see equation (9) 
below — the steady-state population of bound grains will 
be proportionately greater than the unbound population. 
Combining this fact with the tendency for grain size dis- 
tributions to concentrate their collective geometric cross 
section at the smallest sizes, we posit that bound grains 
near the blow-out size, say Sbiow < s < 2sbiow, are re- 
sponsible for much of the observed scattered light. This 
view is consistent with the discussion of particle sizes by 
K05. 

Such grains occupy eccentric orbits, e > 1/3, and are 
disrupted by collisions amongst themselves. The lifetime 
against collisional disruption is 

. . 0 i 1 (AR\ 1/2 

tcol(Sblow ^ S iS 2Sblow) ^ 7 T 7 : ^ 77 I 7 7" ) 

0(a)r LIrtr \ R ) 


(1 - e) 3 / 2 


7 x 10 4 


2/3 

1 — e 


3/2 


yr, 


(9) 


where the effective optical depth r ~ tr(R/ Ai?) 1 / 2 ac- 
counts for the path length ~(i?Ai?) 4 / 2 traversed by a 
grain on a highly elliptical orbit through the birth ring, 
where densities are highest and collisions most likely oc- 
cur. Use of this path length assumes that relative grain 
velocities are of order the local Kepler velocity; this is 
an acceptable approximation for the order unity eccen- 
tricities of interest here. Account of the limited fraction 
of time spent within the torus has also been made, via 
O(a) and (5). 

Compare t co i with the Poynting-Robertson drag time, 
tpR,(sbiow < s < 2 s b iow) = » ' P E(e)R 2 s ~ 2.5 x 10 7 

oIj * 
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2st>lo 


2/3 

1 — e 


1/2 


yr, 


(10) 


which is the time for an orbit of initial periastron R and 
initial eccentricity e to shrink to a point (Wyatt & Whip- 
ple 1950). This is of the same order as the time for a grain 
to have its pericenter be dragged out of the birth ring, for 
A R not much less than R, which is the case here. The 
dimensionless function E(e > 1/3) > 1.9 quantifies the 
decay of orbital eccentricity, and diverges as (1 — e) -1 / 2 . 
Comparison of (9) with (10) shows that as long as e is 
not too close to one — i.e. , for all particle sizes outside 
a tiny interval that just includes Sbiow — grains are re- 
moved from the ring by collisionally cascading down to 
the blow-out size, and PR drag presents only a minor loss 
mechanism. In other words, Fomalhaut’s disk is Type B 
or collision-dominated (Strubbe & Chiang 2006). 

Our estimate of the collisional lifetime t co \ in (9) in- 
forms the duration of our dust particle simulations, in- 
troduced in §3.1.2. 


2.6. Total Dust Mass Versus Total Parent Body Mass 

The mass M d in dust responsible for the scattered light 
is 

M d ~ ^-psr ± RAR ~ 10- 3 M e ) . (11) 

O \ -^blow / 


The mass M p b in the largest parent bodies at the top of 
the collisional cascade is given by the steady-state con- 
dition 


Mpb 

bige 


Md 

tco\ 


(12) 


which implies 


Mpb ~ 3 Mq . 


This is a minimum mass for the disk as a whole because 
still larger bodies may exist which are collisionless over 

t age ■ 

Some workers (e.g., Backman & Paresce 1993) calcu- 
late the mass in parent bodies by explicitly assuming a 
size distribution appropriate for an idealized collisional 
cascade (Dohnanyi 1969) and taking the upper size to 
be some value > 1 km. Not only is it unnecessary to 
specify a size distribution, but it is not justified to adopt 
a specific value for the parent body size without first es- 
tablishing that a typical such body can be collisionally 
disrupted within the finite age of the system. The super- 
kilometer sizes often invoked (e.g., K05) fail this test. 


3. NUMERICAL MODEL 

We devise a dynamical model of the Fomalhaut planet- 
belt system that reproduces approximately some of the 
properties inferred from the HST observations. We com- 
pute the shape of the vertical optical depth profile, tj_ (a), 
of dust particles in the belt and match this profile against 
that of the K05 scattered light model. We seek in partic- 
ular to find those combinations of planet mass and orbit 
that yield an inner edge to the belt of dinner = 133 AU. 

As stressed throughout this paper, the model assumes 
only a single interior planet, an idealization not sup- 
ported by Fom b’s space velocity (§1) or by the Hippar- 
cos data (§4.2). Nevertheless, the planet masses derived 
from the single planet model can usefully be interpreted 
as upper limits (§1). 


The procedure is detailed in §3.1; results are given in 
§3.2; and extensions of the model, including some vali- 
dation tests, are described in §3.3. 

3.1. Procedure 

Our numerical modeling procedure divides into four 
steps, described in the following four subsections, 
§§3.1. 1-3. 1.4. In short, we (1) create a population of 
parent bodies that is stable to gravitational perturba- 
tions from Fom b over f age ; (2) release dust particles from 
stable parent bodies and follow dust trajectories in the 
presence of radiation forces over the collisional lifetime 
t co 1 ; (3) compute the optical depth profile rj_(a) of dust 
particles at the end of t co i, accounting for their size dis- 
tribution; and (4) compare with the K05 scattered light 
model, which serves as our proxy for the observations. 

3.1.1. Step 1: Create Stable Parent Bodies 

Parent bodies (a) execute orbits that are stable to grav- 
itational perturbations over t age , (b) occupy the top of 
the collisional cascade, which by definition implies that 
their orbits are little affected by catastrophic collisions 
for times t < t age , and (c) are large enough that radiation 
forces are negligible. We assume further that (d) the self- 
gravity of the belt is negligible. Thus the problem of sim- 
ulating a realistic parent body is a purely gravitational, 
three-body problem involving the star, planet, and ex- 
terior parent body, where the parent body is treated as 
a test particle. Finding stable test particle orbits is a 
straightforward matter of specifying their initial condi- 
tions, integrating the equations of motion forward for 
~i a ge> and selecting those particles that survive the in- 
tegration. 

Integrations of parent bodies are performed with the 
swift_whm code, written by H. Levison and M. Duncan 
(www. boulder. swri.edu/~hal/swift.html) using the Wis- 
dom & Holman (1991) algorithm. We set the stellar mass 
M* = 2.3 Mq. In each of our five mass models, a planet 
of mass M p i G (0.1, 0.3, 1, 3, 10)Mj resides on a (fixed) 
elliptical orbit of semimajor axis a p i and eccentricity e p i. 
These quantities are listed in Table 1 (placed at the end 
this manuscript); see below for how we relate e p i to a p i, 
and Step 4 (§3.1.4) for how we select a p i given M p \. The 
planet’s orbital plane defines the x-y reference plane for 
the system. The planet’s longitude of periastron w p i = 0. 
At the start of the integration, the planet is located at 
periastron. All orbital elements reported here are stel- 
locentric and osculating unless otherwise noted. 

Each run begins with Nt p = 10 4 test particles and lasts 
10 s yr. Particles are discarded as “unstable” if either 
they approach within a Hill sphere Ah = (/x/3) 1 / 3 a p i of 
the planet, or their distance from the star exceeds 1000 
AU, as a result of gravitational scatterings off the planet. 
Particles that survive until the end of the run are deemed 
“stable” and are used in subsequent steps of our proce- 
dure. The integration timestep is 15000 days, equivalent 
to 5-7% of the planet’s orbital period, depending on the 
model. 

Initial conditions for test particles are as follows. Ini- 
tial semimajor axes a are distributed systematically and 
uniformly between 120 AU and 140 AU (for the 0.1 Mj 
model, the starting particle semimajor axis is 125 AU 
since the planet senrinrajor axis turns out to be a p i = 
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120 AU). Initial eccentricities are set equal to the forced 
values given by the classical, linear, secular theory of 
Laplace-Lagrange (L-L): 


, , _ , *4/2 (“pi/®) 

CyCL) ^forced(^) ^pl 

b 3 M a p i/ a ) 


(13) 


where the b ' s are the usual Laplace coefficients (e.g., Mur- 
ray & Dermott 2000). Initial inclinations i of test par- 
ticles are drawn randomly from a uniform distribution 
between 0 and 0.025 rad. 8 Our maximum inclination 
matches the opening angle of the K05 scattered light 
disk model. Initial longitudes of periastron of all par- 
ticles equal the secularly forced value to = 0, correspond- 
ing to apsidal alignment with the planet; longitudes of 
ascending node O are drawn randomly from a uniform 
distribution between 0 and 27r; and arguments of peri- 
astron to = — f l (so that u = 0). Finally, initial mean 
anomalies M are drawn randomly from a uniform distri- 
bution between 0 and 27r. 

For a given a p i in a given model, the planet’s eccen- 
tricity e p i is such that a test particle at a = 140.7 AU 
acquires a secularly forced eccentricity of e = 0.11, as 
computed using (13). Such parameters are inspired by 
the elliptical orbit fitted by K05 to the brightest portions 
of the belt. The planetary eccentricity so chosen is about 
e p i = 0.12, varying by up to 15% between our five models 
(see Table 1). 

According to L-L, the initial conditions so prescribed 
produce test particle orbits whose eccentricity vectors 
e = (e cos to, e sin to) are purely forced; they have and will 
have no free component (e.g., Murray & Dermott 2000). 
As such, because the planet’s orbit is fixed, belt particle 
orbits also should not vary, at the L-L level of approxi- 
mation. 9 In §3.2.4, we describe the extent to which this 
expectation is borne out. 

Our initial conditions, which comprise nested, apsi- 
dally aligned, purely forced elliptical orbits, are designed 
to reproduce the observed elliptical belt of Fomalhaut 
(Wyatt et al. 1999; Quillen 2006). However, such forced 
orbits are not the only ones that are stable in the vicinity 
of Fom b. In §3.3.3, we experiment with a different set 
of initial conditions that generate another class of stable 
parent body. 


3.1.2. Step 2: Integrate Dust Trajectories 

Having created an ensemble of stable parent bodies, 
we now model the dust generated by such bodies. At 
the end of a 10 8 -yr-long integration from Step 1, we take 
each stable parent body and have it “release” a dust grain 

8 Such an inclination distribution is unphysical because in reality, 
there is zero probability density for finding an object with zero 
inclination. Nevertheless, we adopt our boxcar distribution for 
simplicity. 

9 Laplace-Lagrange truncates the secular disturbing function at 
0 (e 2 , 2 2 ) and so in reality and in numerical integrations, the ec- 
centricities and apsidal angles are still expected to vary somewhat 
with time with our initial conditions, even if the system were purely 
secular. One can avoid such truncation error by employing Gauss’s 
perturbation equations for e and uj and integrating the planetary 
forces over Gaussian wires (e.g., Murray Sz Dermott 2000), thereby 
finding exact forced eccentricities for which e = uj = 0. We skip 
this refinement, in part because the system is not purely secular; 
nearby mean-motion resonances influence the dynamics. 


with the same instantaneous position and velocity as its 
parent’s. Each dust grain’s trajectory is then integrated 
forward under the effects of radiation pressure and PR 
drag. That is, in addition to the gravitational accelera- 
tions from the star and the planet, a dust particle also 
feels a radiative acceleration (see, e.g., Burns et al. 1979) 


GAD f3 

<Uad — ~2 



v r r + v 


(14) 


where r = rr is the vector displacement from the star 
to the grain, v is the velocity of the grain relative to the 
star, and v r = vr accounts for the Doppler shift in stellar 
radiation seen by the grain. We add this radiative accel- 
eration to the Bulirsch-Stoer (B-S) integrator swift_bs, 
written by Levison & Duncan. We prefer to modify the 
B-S integrator over the Wisdom-Holman integrator, since 
the latter is designed to model a dissipationless Hamil- 
tonian system; when PR drag is included, the system is 
dissipative, and it is not obvious to us how we should add 
the radiative perturbations to the symplectic kick-drift- 
kick algorithm. (Nevertheless, adding radiative forces to 
symplectic integrators is standard practice in the liter- 
ature and has been shown to produce accurate results.) 
Though the B-S integrator is slower, it is fast enough 
for our purposes since our integration times for Step 2 
are short (see below). The accuracy parameter “eps” of 
swift_bs is set to 10~ 8 . The modified code was tested 
on test particles having various p’s, producing results for 
radiation blow-out and PR drag in excellent agreement 
with analytic and semi-analytic studies (e.g., Wyatt & 
Whipple 1950). 

From each of the five parent body integrations com- 
pleted in Step 1, we generate eight dust simulations 
in Step 2, each characterized by a single value of /3 € 
(0, 0.00625, 0.0125, 0.025, ..., 0.4). Dust grains released 
with such P’s are mostly still bound to the star, albeit 
on highly eccentric orbits for P approaching 0.4. Bound 
grains contribute substantially, if not predominantly, to 
the scattered light observations: as sketched in §2, be- 
cause tpR, t co \ t esc , the lifetime of bound grains 
is set by destructive interparticle collisions, not by PR 
drag, and the steady-state population of bound grains 
greatly outweighs that of unbound grains, by t co \/t esc . 

Because the dust lifetime is set by collisions, we extract 
dust grain orbits for further analysis in Step 3 after inte- 
grating for a time t co \ . Following our order-of-magnitude 
estimate (9), we set t co i = 10 5 yr. During the integration, 
we discard particles that approach within a Hill sphere of 
the planet or whose distances from the star exceed 10000 
AU (this is a factor of 10 larger than the cut imposed in 
Step 1, because large apastron distances result from the 
onset of radiation pressure and do not necessarily im- 
ply orbital instability). As Table 1 indicates, few if any 
dust particles are discarded in our /3-simulations, with 
the exception of our 10 Mj model. 

In §3.3.2, we test the sensitivity of our results to t co i, 
whose value we know only to within factors of a few. 
In that subsection we also examine whether our results 
change significantly if we model the dust more realisti- 
cally by releasing grains gradually over t co \. 

Figure 1 provides sample snapshots of dust grains and 
their parent bodies projected onto the x-y plane, for our 
lMj model. 
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Fig. 1. — Snapshots of parent bodies (left) and (3 = 0.1 dust 
grains (right), for our lMj model. The cross marks Fomalhaut, 
while the solid circle marks Fomalhaut b. Parent bodies are imaged 
after an integration time of f age + t co l • Dust particles are released 
from parent bodies with zero relative velocity after t age and their 
trajectories integrated forward with f} = 0.1 for t co i. Red ellipses 
correspond to the inner and outer boundaries of the K05 scattered 
light model (a; nner = 133 AU, a ou t er = 158 AU, e = 0.11). Radia- 
tion pressure spreads dust particles outward from where they were 
born, but leaves their inner boundary practically coincident with 
that of parent bodies. 



3.1.3. Step 3: Compute Optical Depth Profile 

If we had a sufficiently large number of dust parti- 
cles, we could simply take a snapshot of each Step 2 
/3-simulation after t co i and count the number of dust par- 
ticles per unit x-y area of the belt. We would thereby 
measure the surface density Ng(x,y ), and by extension 
the vertical optical depth r±(x,y). In practice, we do 
not have enough particles, and such an exercise produces 
noisy results. 

We greatly improve the signal-to-noise by spreading 
each dust particle along its orbit according to how much 
time it spends traversing a given segment of its orbit. 
In other words, we replace each dust particle with its 
equivalent Gaussian wire, and measure the optical depth 
presented by the collection of wires. We refer the Kepler 
elements of a dust particle’s orbit to (1 — /3) times the 
stellar mass; otherwise the elements would not remain 
constant in a two-body approximation. 

First we construct an eccentric grid by partitioning the 
x-y plane into a series of nested, confocal ellipses, all 
having the same eccentricity of 0.11 (K05), and having 
semimajor axes running uniformly from a± = 100 AU to 
on = 220 AU in steps of A a = 0.5 AU. Our goal is to 
measure rj_(ai), the average optical depth between the 
i th ellipse having semimajor axis cq and the {i + l) th el- 
lipse having semimajor axis cq+i = a* + Aa. Each dust 
particle orbit is divided into 1000 segments spaced uni- 
formly in true anomaly from 0 to 27r. Each segment maps 
to a certain location on the grid (i.e., the x-y position of 
each segment falls between two adjacent ellipses on the 
grid). Associated with each segment is an orbital weight, 
equal to the fraction of the orbital period spent travers- 
ing that segment (the sum of all orbital weights for a 
given particle/orbit equals one). The orbital weight for 
each segment is added to its grid location. This process is 
repeated over all segments of all orbits. Finally, at each 
grid location at , the sum of orbital weights is divided 
by the area of the annulus extending from the 2 th ellipse 
to the (2 + l) th ellipse. This yields the relative surface 
density profile Ng{a{), for a simulation characterized by 
(3. 

The various Ng profiles for our 1 Mj model are plotted 
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Fig. 2. — Top panel: Surface density profiles of dust grains for our 
lMj model, computed £ co i = 10 5 yr after release, normalized to the 
peak of the (3 = 0 curve. These profiles are computed by binning 
wire segments on a fixed elliptical grid, as described in the main 
text under Step 3 of our procedure. Profiles shrink vertically and 
widen horizontally with increasing /3, reflecting how increased ra- 
diation pressure spreads particles outward by amplifying apastron 
distances. By contrast, periastron distances are more nearly con- 
served for bound particles, since they always return to the birth 
ring regardless of the strength of radiation pressure. Thus the 
peaks of the Np profiles, which mark the location of the birth ring 
of parent bodies (see left panel of Figure 1), hardly shift with (3. 
Bottom panel: Vertical optical depth profile (solid line) obtained 
by adding together the individual Np profiles (dotted lines), ap- 
propriately weighted according to a Dohnanyi size distribution (see 
equation 15). At a < 160 AU, the profile resembles that of the K05 
scattered light model (dashed line), which itself is an approximate, 
idealized, and non-unique representation of the HST observations. 
Discrepancies at large a > 160 AU are expected, in large part be- 
cause the HST images are too noisy to usefully constrain the K05 
model there. 


in the top panel of Figure 2. The profiles are normal- 
ized to the peak of the N 0 (/3 = 0) profile. Since the 
number of dust particles is practically constant across 
all /3-sinrulations within a given mass model (Table 1), 
the decreasing height of each Ng profile with increas- 
ing [3 simply reflects how dust particle orbits become in- 
creasingly eccentric and elongated with increasing (3 (cf. 
equation 4). In other words, the same number of parti- 
cles is being spread into a disk that extends farther out 
the greater the radiation pressure. At the same time, the 
peaks of the Ng profiles hardly move with increasing (3: 
as long as the grain is still bound to the star, it must 
always return to the same stellocentric distance at which 
it was released, no matter how strongly it feels radiation 
pressure. That release distance is located in the birth 
ring (Strubbe & Chiang 2006) or, more accurately, the 
birth ellipse of parent bodies, near a ~ 130-140 AU; see 
the left-hand panel of Figure 1. 

Once all the Ng profiles are in hand, we combine them 
linearly to produce the total optical depth profile rj_: 
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Tj_ 


y- jy max V0.00625 / 0 \ 9 

j-0 0 max N /3 \ 0.00625 J 


+ N 0 maxNoo ° 625 (l + V2). 

max Nq 


(15) 


The rationale behind this formula is as follows. As Figure 
2 indicates, the maxima of the Np profiles are all situated 
in the birth ring. Following Strubbe & Chiang (2006), 
we posit that the size distribution of grains in the birth 
ring is given by a Dohnanyi (1969) cascade, with differ- 
ential power-law index q = 7/2. For such a power-law 
size distribution, the collective surface area or geometric 
optical depth, evaluated per logarithmic bin in 0, scales 
as 0 q ~ 3 (oc s 3 ~ q ). The two factors multiplying Np in the 
sum in (15) enforce this scaling in the birth ring (whose 
location is traced by the maxima of the surface density 
profiles), and we have adopted the 0 = 0.00625 bin as 
our reference bin. 

The last term proportional to Nq accounts for the opti- 
cal depth contributed by grains having 0 < 0 < 0.00625. 
We take the surface density profile of such grains to be 
given by Nq ; this is a good approximation, as there is 
little difference between Nq. 00625 and N 0 (Figure 2, top 
panel). Since our grid of models for 0 > 0.00625 is log- 
arithmic in 0 — successive bins are separated by factors 
of B = 2 — we extend the same logarithmic grid for 0 < 
0.00625. Then the optical depth coming from all grains 
having 0 < 0.00625, scaled to the optical depth in our 
0 = 0.00625 reference bin, is J/^ 1 (l/J5)^ 9-3 ^ = l + \/2. 


The K05 model is based on fits “by eye.” From the visual 
fits, the uncertainty is about ±1 AU, slightly larger than 
the size of a 2-pixel resolution element (0.1 arcsecond or 
0.77 AU; the images are binned 2x2 before they are 
modeled, to increase signal-to- noise). The error in dinner 
propagates into our constraints on planet mass. 

It is well to appreciate that while our general goal is 
to reproduce the shape of the optical depth profile of 
the K05 scattered light model, that model is itself highly 
idealized, characterized by knife-edge sharp inner and 
outer edges and simple power-law behavior. Fitting by 
eye means the model is at best approximate. And as K05 
caution, only a restricted azimuth of the belt (near their 
quadrant “Q2”) was analyzed in detail to produce their 
fit parameters. Therefore we should neither aim for, nor 
expect, perfect agreement between our dynamical model 
and the K05 model. Our task instead is to use the K05 
model as a guide, to identify robust trends and to rule 
out those regimes of parameter space for Fom b that give 
blatantly poor agreement with observed belt properties. 

Lastly, neither our dynamical model nor the K05 scat- 
tered light model should be trusted at large a > 160 AU. 
At large distance, barely bound grains whose /3’s are ar- 
bitrarily close to the blow-out limit of ~0.5 dominate 
t i_. Our set of 8 /3-simulations lacks the resolution in 
0 to accurately model this outer disk (see Strubbe & 
Chiang 2006 for an analytic treatment appropriate for a 
circular birth ring). Observationally, the HST images at 
stellocentric distances > 158 AU are dominated by the 
sky background; see Figure 3 of K05. 

3.2. Results 


3.1.4. Step 4 : Compare with Observations 

A rigorous comparison with observations would require 
us to produce a scattered light image based on our dy- 
namical model. This is a considerable task. Our tl pro- 
file, combined with the vertical density distribution and 
a grain scattering phase function, would be used to cal- 
culate the direction-dependent emissivity of the belt as 
a function of 3D position. This emissivity model would 
then be ray-traced at a non-zero viewing angle to produce 
a model scattered light image. Various parameters (e.g., 
normalization of t± , grain scattering asymmetry param- 
eter, viewing angle) would need to be adjusted, including 
input parameters from Steps 1 -3 (distribution of initial 
senrinrajor axes, distribution of initial inclinations), to 
produce a good subtraction of the observed image. 

In this first study, we attempt none of this. Instead 
we compare the t± profile given by (15) with the corre- 
sponding optical depth profile of the K05 scattered light 
model, focussing on the one belt property that seems 
most diagnostic of planet mass and orbit: the belt’s in- 
ner edge. The K05 optical depth profile extends from 
dinner = 133 AU to a ou ter = 158 AU and falls as a~ 8 ' 5 
(their fitted volumetric number density of grains falls as 
a~ 9 , while the fitted vertical thickness of the belt in- 
creases as a 0 ' 5 ). In our dynamical modeling, for given 
planet mass M p i, we adjust only the planet semimajor 
axis a p i and repeat Steps 1-3 until the minimum semima- 
jor axis at which t± reaches half its maximum value — the 
“half- maximum radius” — equals ai nner = 133 AU. 

Since our dynamical models are anchored to a^ner, we 
should have a sense of the uncertainty in this parameter. 


In §3.2.1 we give an approximate formula relating the 
planet mass to the width of the chaotic zone, based on 
our five mass models. In §3.2.2, we compare our opti- 
cal depth profiles with that of the K05 scattered light 
model. Based on this comparison we argue against large 
planetary masses for Fom b. In §3.2.3 we argue the same 
point by comparing our model planetary orbits with the 
observed stellocentric distance of Fom b. Finally, the un- 
derlying dynamics of parent bodies and of dust particles 
is discussed briefly in §3.2.4. 

3.2.1. Chaotic Zone Width 

In Figure 3, we overlay the t_l profiles of our five mass 
models together with the optical depth profile of the K05 
scattered light model. As described in §3.1.4, the planet 
senrinrajor axis a p i for each mass model is chosen such 
that the “half-maximum radius” — the smallest senrima- 
jor axis at which t_l attains half its peak value — equals 
K05’s ainner = 133 AU. The a p i’s so determined are listed 
in Table 1 and also annotated on Figure 3. They are such 
that the semimajor axis separation between belt inner 
edge and planet is given by 

ainner &pl = 2.0 /i ^ a p l , (16) 

with less than 6% variation in the coefficient across mass 
models, and where ainner = 133 AU. Because we are 
connecting more directly to the HST observations, our 
Fomalhaut-specifrc coefficient of 2.0 is preferred over the 
smaller coefficients cited by previous works; accuracy 
matters for determining planet mass, whose value scales 
strongly as the 7/2 power of distance measurements. 
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Fig. 3. — Vertical optical depth profiles of our five mass models 
(black lines), overlaid with that of the K05 scattered light model 
(blue dashed line). Parameters for our dynamical models are listed 
in the inset and are such that the “half-maximum radius” — the 
minimum semimajor axis for which r±_ attains half its maximum 
value — equals 133 AU, the innermost semimajor axis of the K05 
model. Models for which M p \ < lMj do equally well in reproduc- 
ing the K05 model. As M p \ increases, the tj_ profiles widen because 
the planet increasingly perturbs dust grains onto eccentric orbits. 
The 10Mj model is probably unacceptably wide. At a > 160 AU, 
neither the dynamical model nor the K05 model is trustworthy; the 
former suffers from poor resolution in (3 , while the latter is limited 
by sky background (see Figure 3 of K05). 


Measured in Hill radii (evaluated using a p i), the sepa- 
ration dinner — dpi ranges from 3.7 to 4.5 in order of 
decreasing M pl . 

3.2.2. Comparison of rj_ Profiles Implies M v \ < 3Mj 

What resemblance there is in Figure 3 between dynam- 
ical and scattered light t± profiles leads us to believe that 
we are on the right track towards understanding the un- 
derlying properties of the Fomalhaut planet-belt system. 
We are especially encouraged when we consider that with 
the exception of a p i, none of our model parameters (e.g., 
grain size index q , distribution of initial semimajor axes) 
has been adjusted from its naive standard value. And as 
we emphasized at the end of §3.1.4, the K05 scattered 
light model is itself highly idealized and approximate, 
and does not represent a unique model of the observa- 
tions. In particular, the K05 model idealizes the inner 
edge as a step function, but smoother profiles also seem 
possible; the degree of smoothness may help to constrain 

q- 

Discrepancies at large a > 160 AU are not serious, since 
both the dynamical and scattered light models are known 
to fail there (§3.1.4). The deficiency in our dynamical 
model can be remedied by having a finer grid in f3 near 
the blow-out value of ~0.5, while improvements in the 
scattered light model await deeper imaging campaigns. 

As Mp\ increases, the t± profiles broaden — see in par- 
ticular the curve for 10Mj. Upon release, dust parti- 
cles find themselves in a weakened stellar potential be- 
cause of radiation pressure. More massive planets more 



Fig. 4. — Snapshots of parent bodies (left) and j3 = 0.1 dust 
grains (right), for our 3Mj model. Compare with the lMj model 
shown in Figure 1. The larger the planet mass, the more dust 
particles are perturbed by the planet into a more spatially extended 
disk. 


readily perturb dust onto more eccentric orbits that ex- 
tend both inside and outside of the birth ring. This is 
further illustrated in Figure 4, which shows a snapshot 
of (3 = 0.1 grains for M v \ = 3Mj. Moreover, as docu- 
mented in Table 1, planetary perturbations are so severe 
for our 10Mj model that about 1/3 of the dust particles 
launched with f3 = 0.05 are discarded as unstable within 
Aoi = 10 5 yr. (Why (3 = 0.05? Larger (3 grains are, upon 
release, repelled immediately away from the planet onto 
highly eccentric trajectories by radiation pressure alone 
and are therefore less likely to be rendered unstable by 
the planet. Smaller (3 grains share essentially the same 
stability properties as the parent bodies.) 

The tj_ profile for M p \ = 10Mj is probably unaccept- 
ably broad: at 140 AU < a < 160 AU, the dynamical 
model predicts too large an optical depth compared to 
the K05 model, by about a factor of two. At these dis- 
tances, a factor of two overluminosity in the belt is not 
easily accommodated, as judged from the error bars on 
the observed surface brightness profile — see Figure 3 of 
K05. This same 10Mj model also produces a tail ex- 
tending inward to a < 120 AU, but the observations, 
whose dynamic range in surface brightness is limited to 
less than 10:1 (see Figure 3 of K05), probably cannot rule 
out such a feature. We have verified that the excessively 
large t± at a > 140 AU follows primarily from the large 
eccentricities acquired by [3 « 0.2-0. 4 dust particles from 
planetary perturbations. By contrast, the M v \ < lMj 
models, which produce practically identical t_l profiles, 
appear compatible with the K05 model, given the various 
limitations of the latter. 

The low-mass models having M p \ < lMj have inner 
edges that are practically identically sharp. A mea- 
sure of the sharpness is the distance over which t± rises 
from the half-maximum radius to the semimajor axis at 
which r l peaks, normalized by the half-maximum radius. 
This fractional distance is S = 4.5AU/133AU = 0.034. 
Quillen (2006) recognizes that in fitting the observed 
drop-off in surface brightness for a belt that is inclined to 
our line of sight, a trade-off exists between the belt’s ver- 
tical and radial density profiles. Either the belt can have 
a finite vertical thickness and be knife-edge sharp in the 
radial direction (as in the K05 scattered light model); or 
have zero vertical thickness and drop off gradually in the 
radial direction; or be characterized by some intermedi- 
ate combination. In this context, our measure for the 
fractional radial drop-off distance, <5 = 0.034, compares 
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promisingly with the fractional vertical drop-off distance 
inferred by K05, 2 H/R = 0.025. In the future, we will 
have to adjust both radial and vertical drop-off distances 
to better reproduce the scattered light observations. The 
vertical thickness of the belt appears to be fairly easily 
adjusted in our model, given that the inclinations of our 
dust particles are mostly unchanged from the assumed 
initial inclinations of our parent bodies (see Figure 7 and 
related discussion in §3.2.4). 

In the bottom panel of Figure 2, we plot as dotted 
curves the separate terms in equation (15) that add up 
to tj_. The terms corresponding to larger (3 (or smaller 
grain size s oc 1//3) dominate, as a consequence of our 
assumption that the grain size distribution in the birth 
ring follows a Dohnanyi law. That law apportions the 
bulk of the geometric surface area in the smallest grains. 

3.2.3. Fom b ’s Current Stellocentric Distance Implies 

M p i < 3Mj 

Each of our mass models specifies a certain orbit for 
Fom b (Table 1) that is tuned, through multiple itera- 
tions of Steps 1-4, to generate the observed ellipticity 
of the belt and to yield a half-maximum radius equal to 
K05’s dinner = 133 AU. If the apocentric distance of a 
model orbit is less than the observed stellocentric dis- 
tance of Fom b, then that model can be ruled out. 

Currently, only the distance between Fom b and its 
host star projected onto the sky plane is known. In 2006, 
that projected distance was 97.6 AU. Inferring the true 
stellocentric distance requires that we de-project the or- 
bit of Fom b. But with only two epochs of observation, 
a unique de-projection is not possible. Nevertheless, we 
can perform a rough de-projection by making a few as- 
sumptions. We assume that the planet’s orbit is oriented 
such that its senrinrajor axis lies in the plane of the sky 
and is parallel to the line joining the observed belt ansae. 
We also assume that the inclination of the planet’s orbit 
to the sky plane equals 65.6°, the same as that inferred 
for the belt orbital plane (K05). These assumptions are 
reasonable insofar as they are self-consistent with the sin- 
gle planet hypothesis, that Fom b dominates the belt’s 
dynamics. The single planet hypothesis predicts that 
belt and planet orbits are apsidally aligned and copla- 
nar. 

The resultant de-projected stellocentric distance of 
Fom b in 2006 is 119 AU, with a systematic error of 
probably no more than a couple AUs. Such a distance 
argues against our 10Mj model, for which the apocen- 
tric distance of Fom b is 107 AU. The discrepancy is 
depicted in Figure 5. That same figure shows that the 
3Afj model is also inconsistent, but only marginally so, 
and without a proper de-projection, we hesitate to rule 
it out. Lower masses M p \ < lMj are, by contrast, easily 
compatible. These findings reinforce those based purely 
on a comparison of the optical depth profiles (§3.2.2). 

3.2.4. Parent Body and Dust Particle Dynamics 

In Figure 5, we supply sample histograms of time- 
averaged semimajor axes of stable parent bodies gener- 
ated in Step 1. The time average is performed over a 10 5 - 
yr-long window (with (3 = 0) following each 10 8 -yr-long 
integration. Intriguingly, parent bodies appear evacu- 
ated from exterior mean-motion resonances established 
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Fig. 5. — Histogram of time-averaged semimajor axes of parent 
bodies that survive for 10® yr in the vicinity of Fom b. The bin 
width is 0.1 AU, and the time average is performed over 10 s yr. 
Each panel corresponds to a different mass model, as annotated. 
The black circle in each panel marks the semimajor axis of Fom 
b, with error bars extending from the model orbit’s pericentric 
distance to its apocentric distance. Only model orbits correspond- 
ing to M p i < ltfj are consistent with the estimated de-projected 
stellocentric distance of Fom b in 2006 (dotted vertical line). Sta- 
ble parent bodies are located outside the planet’s chaotic zone, at 
semimajor axes greater than the planet’s by about 2// 2 ' ' a p i . In- 
side the chaotic zone, first-order mean-motion resonances overlap 
and particle orbits are short-lived. Outside the chaotic zone, par- 
ent bodies reside stably on secularly forced eccentric orbits, with 
occasional gaps located at mean-motion resonances as indicated. 
The gaps are evacuated for reasons likely analogous to why the 
solar system’s Kirkwood gaps are empty of asteroids. 


by Fom b, even outside the planet’s main chaotic zone. 
The reasons for this are likely analogous to why the so- 
lar system’s Kirkwood gaps, located at Jupiter’s interior 
mean-motion resonances, are empty of asteroids (Wis- 
dom 1982, 1983, 1985). Study of this phenomenon, which 
depends on the non-zero eccentricity of the planet’s orbit, 
is deferred to future work. 

Figure 6 plots the eccentricity vectors of the parent 
bodies at the end of 10 s yr (black points). While they 
remain clustered near their initially purely forced values, 
there is a dispersion that is not predicted by L-L: the par- 
ent bodies acquire free eccentricities despite having none 
to start with. The more massive the planet, the greater 
the free eccentricities that develop. This same behavior 
was found by Quillen (2006) and Quillen & Faber (2006), 
who attributed it to forcing by a mean-motion resonance 
just outside the chaotic zone. 

Figure 6 also describes how dust particles (colored 
points) born from parent bodies acquire free eccentric- 
ities. Initial free eccentricity vectors are distributed 
roughly axisymmetrically about the forced eccentricities. 
The initial free phase depends on the orbital phase of re- 
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Fig. 6. — Locations in complex eccentricity space of all parent 
bodies at t = 0 (green bar) and those parent bodies that survive for 
t = 10® yr (black points). The locus of survivors is not exactly cen- 
tered on the locus of initial conditions because the survivors are all 
at large semimajor axes (smaller forced eccentricities), and because 
of inaccuracies in the L-L theory which was used to determine the 
initial conditions. Surviving parent bodies remain approximately 
apsidally aligned with the planet’s orbit, deviating by less than 
±15° in most cases. The deviations, i.e., the dispersions in free 
eccentricity, increase with M p \. This same behavior was found by 
Quillen (2006). However, contrary to that work, we find that the 
increased dispersion does not necessarily imply a more spatially 
diffuse inner edge to the belt; see Figure 8 and §3.3.1. Also shown 
are 0 = 0.1 dust particles 10 4 yr after release (red points), and 
those same dust particles 10 s yr after release (blue points). The 
trajectory of a typical dust particle is shown sampled every 10 3 yr, 
starting from release. Note the large eccentricity variations for the 
3 Mj model. 


lease. For example, a dust particle released at the parent 
body’s periastron gains a free eccentricity vector in the 
+k direction (total eccentricity increases), while release 
at apastron yields a free eccentricity in the — k direction 
(total eccentricity decreases, at least for (3 not too large). 
Though these radiation-induced free eccentricities of dust 
grains are much larger than the free eccentricities ac- 
quired by parent bodies, they do not necessarily lead to 
increased blurring of the belt inner edge, because the 
semimajor axes of dust grains shift correspondingly out- 
ward by radiation pressure as well (cf. equation 5). We 
have verified that the various Np surface density profiles 
for M p i < lMj have comparably sharp inner edges for 
all /3, as gauged using our fractional width parameter 5. 
For further discussion of what influences the sharpness 
of the inner edge, see §3.3.1. 

Finally, in Figure 7 we examine how the orbital in- 
clinations of stable parent bodies change after 10 8 yr. 
Laplace-Lagrange predicts that the mutual inclination 
between planet and particle is conserved, and indeed 
most inclinations are little altered from their initial val- 
ues (which span up to 0.025 rad). Fewer than 10% of all 
surviving parent bodies have final inclinations that ex- 
ceed initial inclinations by more than 0.0025 rad. Excited 


Fig. 7. — Changes in orbital inclinations of surviving parent bod- 
ies, evaluated after 10 s yr. Solid circles denotes the planet. Only 
small percentages of stable parent bodies have their inclinations in- 
creased by the planet by more than 2.5 X 10“ 3 rad (0.14 deg); these 
percentages are indicated in each panel. For comparison, the initial 
distribution of inclinations extends uniformly from 0 to 2.5 X 10" 2 
rad. Those few objects that have their inclinations amplified are 
localized to mean-motion resonances, labelled only for the middle 
panel. 


bodies are localized to mean-motion resonances. Even in 
the vicinity of a resonance, the fraction of bodies that 
have their inclinations pumped is small. For example, 
in our lMj model, only 10% of all stable parent bodies 
having time-averaged semimajor axes between 132 and 
133.1 AU — particles right at the inner edge of the belt, 
in the vicinity of the 4:3 resonance — have final inclina- 
tions exceeding initial inclinations by more than 0.0025 
rad. The near constancy of inclination should simplify 
future modelling efforts: whatever vertical thickness of 
the belt is desired to match the scattered light images 
can be input into the dynamical model as an initial con- 
dition, in Step 1. 

3.3. Extensions and Refinements 

In §3.3.1 we reconcile our results on the sharpness of 
the inner belt edge with those of Q06. In §3.3.2 we ex- 
amine how robust our optical depth profiles are to un- 
certainties in t co i and to our simplifying assumption that 
the profile is adequately simulated by releasing grains 
from parent bodies at a single time. In §3.3.3 we ex- 
periment with different test particle initial conditions to 
see whether they might yield superior fits to the obser- 
vations. 

3.3.1. Sharpness of Inner Belt Edge 

In Figure 3, we found that the sharpness of the belt’s 
inner edge, as gauged by our measure S (see §3.2.2), 
hardly varied with M p \ < lMj. This finding is seemingly 
at odds with that of Q06, who gauges sharpness using the 
velocity dispersion of parent bodies at the boundary of 
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the chaotic zone, and who finds that it increases smoothly 
as /x 3 / 7 . As a consequence of this relation, Q06 concludes 
that M p i cannot exceed ~7 x 10 -5 M* = 0.2Mj and still 
have the belt edge be as sharp as the observations imply. 
Indeed we also found in Figure 6 that the velocity disper- 
sion of parent bodies, as indicated by their spread in free 
eccentricities, increased smoothly with M p i, in apparent 
agreement with Q06. Yet the smoothly growing velocity 
dispersion is not reflected in the relative sharpness of our 
optical depth profiles across mass models. 

How can this be? It might be thought that the dis- 
crepancy arises because Q06’s calculation of the velocity 
dispersion pertains to parent bodies, while our calcula- 
tion of t± involves dust. While as a point of principle our 
calculation would be preferred because the observations 
are of dust and not of parent bodies, this explanation 
does not get at the heart of the problem, as we find the 
same invariant sharpness with planet mass characterizing 
the surface density profiles of our parent bodies. 

The answer instead is that the sharpness of the in- 
ner edge does not depend only on the velocity dispersion 
of particles located strictly at the chaotic zone bound- 
ary. Particles located at some radial distance from the 
boundary, further interior to the belt, also contribute 
to the sharpness. That is because sharpness is a rela- 
tive quantity, measured relative to the maximum of t±, 
and this maximum does not occur exactly at the chaotic 
zone boundary. Sharpness is appropriately measured as 
a relative quantity, since the observations are limited in 
dynamic range: as Figure 3 of K05 indicates, only the 
maximum in surface brightness, and values greater than 
~10% of the maximum, are measurable. 

To illustrate our point, we plot in Figure 8 the sur- 
face densities of stable parent bodies for two mass mod- 
els, 0.3Afj and 10Mj. For the comparison with Q06 to 
be fair, we must analyze only the parent bodies, since 
the conclusions of Q06 regarding inner edge sharpness 
pertain to collisionless, radiation-free particles; in other 
words, the test particles of Q06 are our parent bodies. 
From Q06 we would expect that the 10Mj model pro- 
duces an inner edge that is (10/0. 3) 3 / 7 = 4.5x more 
diffuse, but the bottom panel of Figure 8 shows that 
this is not the case; in fact, if anything, the 0.3Mj pro- 
file appears more diffuse. In the top panel of Figure 8 
we show the same two models except that we include 
only particles having the smallest stable semimajor axes: 
a < 132 AU for M pl = 0.3Mj and a < 136 AU for 
M p i = 10Mj. These profiles, which more strictly sam- 
ple the chaotic zone boundary, are more consistent with 
Q06 — the 0.3Mj profile is steeper than the 10Mj profile. 
We conclude that sharpness cannot be reliably computed 
without including contributions from particles that are 
located some distance from the edge. Sharpness cannot 
be calculated as if it were a local quantity specific to the 
minimum stable semimajor axis; the shape of the ob- 
served surface brightness profile reflects an amalgam of 
semimajor axes. 

3.3.2. Variations in t co \ and Gradual Release of Grains 

According to our standard procedure, grain surface 
densities and optical depths from our /3-simulations are 
extracted after an integration time of t co i = 10 5 yr, a 
value inspired from our order-of-magnitude estimate (9) 
of the collisional lifetime. In the top and middle pan- 



Fig. 8. — Reconciling our findings on the sharpness of the belt 
inner edge with those of Q06. Top panel : Surface density profiles 
of stable parent bodies for two mass models, 0.3Mj and 10Mj, 
including only those bodies having the smallest semimajor axes as 
indicated (cf. Figure 5). Only parent bodies are considered to 
compare fairly with Q06, whose conclusions regarding inner edge 
sharpness pertain to collisionless, radiation-free test particles. In 
qualitative agreement with Q06, the 10Mj model yields a more 
diffuse inner edge, a consequence of that model’s larger velocity 
dispersion at the chaotic zone boundary (cf. Figure 6). Bottom 
panel : Surface density profiles of the same two models with no 
restriction on semimajor axes. Now the inner edges are of similar 
sharpness — in fact the 10Mj profile appears slightly sharper than 
the 0.3Mj model — indicating that sharpness is not uniquely related 
to edge velocity dispersion, contrary to the implicit assumption of 
Q06. 


els of Figure 9, we demonstrate that our results are not 
sensitive to uncertainties in f co i. We present test results 
for (3 = 0.2 and (3 = 0.4 since those cases contribute 
most to t l for our assumed Dohnanyi size distribution. 
Varying t co \ from 3 x 10 4 yr to 3 x 10 5 yr produces prac- 
tically identical results for the surface density of grains. 
Only for the f3 = 0.2, t co[ = 3 x 10 5 yr run is there a 
slight ~1 AU shift of the surface density profile, a con- 
sequence of Poynting-Roberton drag. Such drag affects 
(3 < 0.2 grains even less. The highest /? grains are also 
relatively immune from pericenter decay — see the middle 
panel of Figure 9 for the case (3 = 0.4- because of the 
large eccentricities and semimajor axes induced by radi- 
ation pressure upon release (Wyatt & Whipple 1950, see 
also our equation 10). 

We can also check whether our simple procedure of 
releasing grains at a single time is sound. In reality, 
the belt at any given moment will contain grains having 
a variety of ages, ranging from 0 (just released grains) 
to tco\ (grains just about to be shattered to sizes small 
enough to be blown out by radiation pressure). We better 
simulate the gradual release of grains by adding together 
surface density profiles computed from grains released 
at multiple times. We divide the interval f C oi into 10 
release times, uniformly spaced by A t = t co i/10, and 
generate separate integrations for each time. That is, for 
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Fig. 9. — Top panel : Surface density profiles for (3 = 0.2 particles 
in our lMj model, as a function of t co \, whose value we only know 
to order of magnitude. Aside from a small inward shift caused by 
PR drag for t co \ = 3 x 10 5 yr, the profiles are not sensitive to our 
uncertainty in t co \. Middle panel : Same as top panel, except for 
(3 = 0.4. Bottom panel: Testing our simplifying assumption of a 
single release time for dust grains. Allowing grains to be released at 
ten uniformly spaced times between t = 0 and t = 10 5 yr produces 
no discernible difference from our standard model, which assumes 
a single release time of t = 0. 


the i th release time, we integrate the stable parent bodies 
(generated from Step 1) forward with f3 = 0 for iAt , and 
then continue integrating the released dust grains with 
/ 3 ^ 0 for (10 — i) At. The resulting superposition of 
integrations, for ft = 0.2 and f co i = 10 5 yr, is shown in the 
bottom panel of Figure 9; it is nearly indistinguishable 
from our standard result. We have checked that ring 
shape is independent of grain age t for I/O, t < t co \, 
because f co i is too short a time for grains to evolve away 
from their birth orbits (say by PR drag), and 1/12 is the 
timescale for grains to phase mix. 

The experiments described in all three panels of Fig- 
ure 9 also show that our standard surface density and 
optical depth profiles do not reflect dust features that 
vary with the orbital phase of the planet. In varying f co i 
and superposing snapshots taken at different times, we 
are extracting surface densities corresponding to random 
orbital phases of the planet. 

3.3.3. Resonant Particles as Another Class of Stable 
Parent Body 

All our parent bodies are initialized with purely 
forced eccentric orbits, as calculated using the Laplace- 
Lagrange secular theory. Their orbits after 10 8 yr resem- 
ble their initial ones, with the addition of a small free 
component. 

Here we try a different set of initial conditions: os- 
culating e = 0. According to L-L, this corresponds to 
assigning particles free eccentricities equal in magnitude 
to their forced eccentricities (see, e.g., Murray & Der- 
rnott 2000). The top panel of Figure 10 documents the 
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Fig. 10. — Experimenting with initially zero eccentricities for 
parent bodies, for the case M p \ = 0.3Mj. Top panel: Histogram 
of semimajor axes of parent bodies that survive for 10 8 yr, time- 
averaged over 10 5 yr. In stark contrast to our standard model 
(Figure 5), survivors only inhabit mean-motion resonances, which 
afford them protection from the close planetary encounters that 
would otherwise result from the large eccentricities that develop 
secularly. Bottom panel: Optical depth profile of dust generated 
by resonant parent bodies. It is far too broad to agree with the K05 
model. We conclude that while resonant parent bodies can exist 
in principle, their population in reality must be small compared to 
that of parent bodies on nearly purely secularly forced orbits. 


resultant time-averaged semimajor axes of stable par- 
ent bodies (those that survive for 10 8 yr), for the case 
M p i = 0.3Mj. All other initial conditions apart from the 
test particles’ eccentricities are specified the same way as 
for our standard 0.3Mj model. Remarkably, comparing 
the top panels of Figures 10 and 5, we find the distribu- 
tions of stable semimajor axes are nearly complementary. 
Instead of being cleared out of mean-motion resonances, 
stable parent bodies inhabit them exclusively when ini- 
tial eccentricities equal zero. 

The resonances — which include the 7:6, 6:5, 5:4, 9:7, 
and 4:3, and which are of eccentricity-type — protect the 
particles from close encounters with the planet. Quali- 
tatively, the eccentricities and apsidal angles behave as 
L-L predicts: while e cycles from 0 through max(e) = 
2ef 0 rced ~ 0-22 back to 0, w regresses from n/2 to 
— 7r/2 (the evolution of uj is discontinuous since e passes 
through 0). The large maximum eccentricities, which re- 
sult because the free eccentricities are of the same magni- 
tude as the forced eccentricities, put particles in danger of 
close planetary encounters, especially when e = max(e) 
and ui = 0. Under these conditions, for a semimajor axis 
of, say, a = 128 AU, the particle’s pericenter encroaches 
within ~2 AU « 0.5i?n of the planet’s pericenter. 

However, thanks to resonance, conjunctions occur only 
at special orbital phases and close encounters do not oc- 
cur. For the circumstances just described, during the 
phase that the particle attains maximum eccentricity, we 
have observed in our numerical integrations that conjunc- 
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Fig. 11. — Same as Figure 1, except for resonant parent bod- 
ies. The left panel shows how resonant parent bodies avoid close 
encounters with the planet. The entire pattern corotates with the 
planet. Dust grains released from these resonant parent bodies 
have large eccentricities and yield an optical depth profile too broad 
to match that of K05; see Figure 10. 


tions never occur at periastron. Because of the 7:6 reso- 
nance, they occur instead 180° degrees away, at apoapse, 
when the bodies are well separated by about 27 AU. Fig- 
ure 11 provides a snapshot of stable resonant parent bod- 
ies, showing how they avoid approaching the planet. 

Can such resonant parent bodies be present in Fomal- 
haut’s belt? Not in significant numbers compared to the 
non-resonant population. The resonant bodies develop 
more eccentric orbits, and consequently the dust they 
produce is more spatially extended. From the bottom 
panel of Figure 10 (see also the right-hand panel of Fig- 
ure 11), it is clear that the optical depth profile of dust 
released from resonant parent bodies is far too broad to 
match that of K05. (Our procedure of calculating opti- 
cal depths by smoothing particles over their orbits is not 
correct for resonant objects, since the smoothing ignores 
their special orbital phase relationships with the planet. 
However the error accrued is small, since t± is domi- 
nated by dust particles having (3 > 0.1. Such dust par- 
ticles, upon release, have their semimajor axes increased 
by > 10% by radiation pressure, and are thus removed 
from the resonances inhabited by their parents.) 

In §4.3 we discuss how the parent bodies might have 
come to occupy nearly purely secularly forced orbits and 
to avoid the resonant orbits. 

4. SUMMARY AND DISCUSSION 

We review our main results in §4.1; sketch the effects 
of other planets apart from Fom b in §4.2; and discuss 
possible origins of Fom b and the belt in §4.3. 

4.1. Summary 

Fomalhaut b is the first extrasolar planet candidate to 
be directly imaged at visible wavelengths and to have its 
orbital motion around its host star measured. Surprises 
have been immediate: Fomalhaut b has an unusually 
large orbital radius of more than 110 AU (cf. Lafreniere 
et al. 2008), and a visual (0.45-0.7 /im) luminosity that 
is not only 1 2 orders of magnitude greater than atmo- 
spheric models anticipated, but also time variable. Nev- 
ertheless, a chain of arguments based on comparing the 
observed photometry with model exoplanet atmospheres 
leads Kalas et al. (2008, K08) to infer that the mass of 
Fomalhaut b must be less than about 3Mj. 

The Fomalhaut system is all the more remarkable for 
offering a means independent of model spectra to get at 
Fom b’s mass: the star is encircled by a belt of dust whose 


geometry is, in principle, sensitive to the mass and orbit 
of Fom b. At a system age of ~200 Myr, detritus from 
the formation of the Fomalhaut planetary system still re- 
mains. If Fom b is the sole sculptor of this debris — but 
see §1 and §4.2 for reasons that it might not be — then 
the observed intrinsic ellipticity of the ring would owe its 
origin to secular forcing by Fom b, which itself would re- 
side on an apsidally aligned and similarly eccentric orbit 
to the belt’s (Wyatt et al. 1999; Quillen 2006). Under 
the single planet assumption, another feature of the belt 
influenced by Fom b would be its inner edge. The ob- 
served sharpness with which the belt truncates would re- 
flect the sharp divide between stability and chaos at the 
boundary of the planet’s chaotic zone, inside of which 
first-order mean-motion resonances overlap and particle 
orbits are short-lived (Wisdom 1980; Quillen 2006). Par- 
ticles inside the zone quickly evolve onto planet-crossing 
orbits and thereafter are perturbed onto escape trajec- 
tories. 10 Identifying the belt inner edge with the chaotic 
zone boundary relates the distance between the planet 
and the belt edge to the planet-star mass ratio. 

Based on these and other theoretical ideas, we have 
built a realistic dynamical model of the Fomalhaut 
planet-belt system, under the single planet assumption. 
Our goal is to calculate the spatial distribution of dust 
generated from the collisional comminution of larger par- 
ent bodies, and to compare our dust maps with the HST 
scattered light observations. The model begins by using 
numerical integrations to establish an annulus of par- 
ent bodies — a.k.a., the “birth ring” (Strubbe & Chiang 
2006) — that is stable for 100 Myr against perturbations 
by Fom b. Dust particles are released from these dy- 
namically stable parents, and their trajectories followed, 
taking additional account of stellar radiation forces, for 
a dust collisional lifetime of t co \ = 0.1 Myr. The orbits 
of dust particles at the end of a t co i-long integration are 
assumed to represent well those of actual dust particles, 
which in reality have a range of ages extending up to 
~f c oi- We have subjected this assumption to a few tests 
and found it to hold. Final dust particle orbits, com- 
puted for a range of radiation /3’s (force ratio between 
stellar radiation pressure and stellar gravity), are con- 
verted into orbit-averaged maps of relative surface den- 
sity (number of grains per unit face-on area of the belt). 
A Dohnanyi (1969) grain size distribution, appropriate 
for a quasi-steady collisional cascade, is assumed to hold 
in the birth ring, where the dust surface density is high- 
est and collision rate is greatest. This assumption deter- 
mines how we weight and add the surface density maps to 
produce a map of (relative) vertical optical depth. That 
optical depth map — or rather its azimuthally averaged 
version, the variation of optical depth with semimajor 
axis- is compared with the optical depth profile of the 
Kalas et al. (2005) scattered light model, which itself 
represents an idealized and approximate fit to the HST 
images. 

A conservative result of our dynamical model is that 
M p i < 3Mj. This conclusion, that Fom b must be 
of planetary mass, is entirely independent of Fom b’s 
photometry — and its uncertain interpretation. Our re- 

10 We estimate that planet crossing takes ~10 5 (10 — 3 //r) 4 U yr 
and ejection takes ~10 7 (lCU 3 //r) 2 yr, for particles halfway between 
the planet and the edge of the chaotic zone. 
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suit stems from two simple and robust trends, neither of 
which involve inner edge sharpness, in contrast to Q06. 
First, as M p \ increases, dust particles are increasingly 
perturbed by the planet onto more eccentric orbits, ren- 
dering the dynamical profiles too broad compared to the 
scattered light profile. A mass of 10Mj yields a belt that 
is about twice as bright between 150 and 160 AU as the 
observations allow, while lower masses M p \ < 3Mj give 
optical depth profiles that we feel agree adequately well 
with the Kalas et al. (2005) scattered light model. Sec- 
ond, given the observed location of the inner edge of the 
belt, larger mass planets have necessarily smaller orbits 
located farther interior to the belt, and smaller orbits 
may be incompatible with the observed stellocentric dis- 
tance of Fom b. A mass of 10Mj requires an orbit whose 
apocentric distance is 107 AU, falling well short of our 
estimated de-projected distance for Fom b of 119 AU. 
By contrast, for M p \ < lMj, the model apocentric dis- 
tances are > 122 AU. While 3Mj yields an orbit whose 
apocenter lies at 115 AU and is nominally incompati- 
ble, uncertainties in the observed de-projected distance, 
probably amounting to a few AU, preclude us from ruling 
out this mass. Thus, erring on the safe side, we conclude 
that M p i < 3Mj. 

Our findings agree in broad outline with those of 
Quillen (2006), whose prediction that the belt inner edge 
be located at the boundary of the planet’s chaotic zone 
appears vindicated by the discovery of Fom b (modulo 
the nominal finding that belt and planet orbits are apsi- 
dally misaligned —whether the misalignment is real will 
have to await multi-epoch astrometry of Fom b). We 
diverge from Quillen (2006) in how we determine Fom 
b’s mass. Among the improvements we make are: (1) 
we draw a clear distinction between unobservable par- 
ent bodies and observable dust grains, and rely on the 
latter when comparing with the HST scattered light ob- 
servations; (2) we include the effect of stellar radiation 
pressure, significant for dust grains; (3) parent bodies 
are screened for dynamical stability over the age of the 
system; and (4) grain-grain collisions are recognized as 
destructive, and therefore the duration of each of our 
dust particle integrations is necessarily limited by the 
collision time. Our upper mass limit of 3Mj follows, in 
part, from comparing theoretical and observed optical 
depth profiles of dust, computed globally over all space, 
and from noting how steep those profiles are both inside 
and outside of the location of peak optical depth. By 
contrast, the upper mass limit derived by Quillen (2006), 
~0.2Mj, follows from an analysis of the radiation-free dy- 
namics of collisionless particles — essentially, the dynam- 
ics of parent bodies — local to the chaotic zone boundary, 
and from the assumption that the local velocity disper- 
sion of such bodies determines the sharpness of the inner 
belt edge. Our calculation of the upper mass limit is 
preferred because the HST observations are of dust, not 
of parent bodies, and also because we have shown that 
edge sharpness is better computed using a global model 
such as ours. 

Three Jupiter masses for the mass of Fom b is an upper 
limit in still another sense: our results are based on the 
assumption that Fom b alone sculpts the belt. Roughly 
speaking, the more planets that are present, the more 
orbits are chaotic, and the more effectively small bodies 
and dust are gravitationally scoured. Thus the observed 


inner edge of Fomalhaut’s belt is compatible with a mass 
for Fom b lower than that computed under the single 
planet assumption. It is heartening to see that our pre- 
ferred mass range of M p \ < 3Mj supports that inferred 
from the spectral models. 

Perhaps the biggest deficiency of our model lies in our 
crude treatment of grain collisions. In setting a dust 
particle’s initial position and velocity equal to that of its 
parent body, we ignore collisional dissipation and redi- 
rection of orbital kinetic energy. We also neglect the fact 
that grinding the largest parent bodies down to dust re- 
quires multiple collisions, and that radiation effects can 
start manifesting in the middle of the collisional cascade. 
For example, a particle for which radiation effects are sig- 
nificant, say which has f3 = 0.4, can be born from a par- 
ent body for which radiation effects were also significant, 
say which had [3 = 0.2. Still, our simple prescription of 
releasing dust grains having /3 > 0 from bodies having 
/3 = 0 is not without justification. Collisional energy dis- 
sipation should, on average, dampen free eccentricities 
but not forced eccentricities (see §4.3); thus the mean 
elliptical shape of the belt is expected to be preserved. 
Grain ejection velocities relative to the parent are dis- 
tributed isotropically and should therefore not bias our 
results, though they will produce larger free eccentrici- 
ties than our model predicts. Larger free eccentricities 
will result in smoother optical depth profiles: a blurring 
of the belt (see §4.3 for further discussion of free eccen- 
tricities). Finally, radiation effects are predominantly felt 
over only the last decade in grain size above the blow-out 
size, i.e. , only after the penultimate collision just prior to 
the final collision resulting in blow-out, for collisions that 
are strongly disruptive. The next generation of models 
should test these assertions, in addition to considering 
qualitatively different physics (e.g., Yarkovsky drag, and 
gas-particle interactions. 11 ) 

4.2. Other Planets In Addition to Fom b 

As discussed in §1, the observed space velocity of Fom 
b is nominally inconsistent with its orbit being apsidally 
aligned with that of the belt, contradicting a basic pre- 
diction of the hypothesis that Fom b is solely responsible 
for the morphology of Fomalhaut’s debris disk. The final 
word on orbits must await more epochs of astrometry 
and more realistic assessments of systematic uncertain- 
ties (in, e.g., frame registration between epochs). Even if 
the apsidal misalignment proves real, and even if follow- 
up observations constrain Fom b’s mass to be less than 
that required to make a significant contribution to shap- 
ing the belt, ideas of chaotic zone clearing and secular 
forcing of eccentricity may still be relevant for the Foma- 
lhaut system, provided there are perturbers in addition 
to Fom b. 

We have already argued that multiple planets are com- 
patible with the inner edge to Fomalhaut’s debris disk; 
more planets help eject more material from inside the 
hole. Multiple planets are also compatible with the ob- 
served belt eccentricity. In Laplace-Lagrange theory, 
the observed forced eccentricity vector of a belt parti- 
cle equals the vector sum of n eccentricity vectors forced 

11 Fomalhaut’s closest analogue may be AU Mic, insofar as both 
have birth ring morphologies and similar optical depths (Strubbe 
& Chiang 2006). Gas has not been detected in AU Mic (Brandeker 
et al., 2008 Spitzer Science Conference Poster #81). 
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by n planets. Because the individual eccentricity vec- 
tors precess at frequencies that depend only on the fixed 
masses and semimajor axes of the planets (these are the 
fixed eigenfrequencies of the linear theory), and because 
belt particles having similar semimajor axes have similar 
vector decompositions, we expect the forced eccentricity 
vectors of belt members to remain similar to one another 
over time. Thus a narrow belt can maintain a global 
mean eccentricity in the presence of multiple planets, 
though that mean eccentricity will oscillate with time. 
All these considerations can be accommodated as neces- 
sary into our modelling procedure. 

If future astrometry confirms a significant apsidal mis- 
alignment between planet and belt, then at least one 
other, as yet unseen planet would be implicated. In that 
case, because the belt’s forced eccentricity is a vector 
sum, Fom b’s eccentricity could either be lower or higher 
than the eccentricity we have calculated, depending on 
the apsidal orientation(s) of the other planetary orbit(s). 

Intriguingly, over its three-year mission, the Hippar- 
cos satellite observed Fomalhaut to have an “anomalous” 
proper acceleration of 6.6 milliarcsec/yr , of marginal 
(about 2a) significance. If real, such a quasi-steady accel- 
eration might be caused by a companion whose orbital 
period is longer than ~3 years, or equivalently whose 
stellocentric distance r > 3 AU. Equating the measured 
acceleration with ( GM / r 2 )/d , where M is the perturbing 
mass and d = 7.7 pc is the distance to Fomalhaut, we see 
that Fomalhaut might also be harboring a ~30Mj brown 
dwarf at a distance r ~ 5 AU. (Larger r implies larger M 
and such solutions are probably ruled out by observation. 
Clearly Fom b cannot be responsible for the Hipparcos 
acceleration.) Compared against the influence of Fom b 
as we have computed it in this paper, such a brown dwarf 
would contribute more than 10% to the forced eccentric- 
ity of a belt particle, if the brown dwarf’s eccentricity 
> 0 . 2 . 

The above lines of evidence for additional planets per- 
turbing the belt are tenuous. Given the observed prox- 
imity of Fom b to the belt, the opposing case can be 
made that Fom b dominates the forced eccentricities of 
belt particles. A possible analogy would be Neptune and 
the Kuiper belt. Despite the existence of as many as 
four giant planets in our solar system, the forced inclina- 
tion and eccentricity vectors of Kuiper belt objects are 
largely determined by the nearest planet, Neptune (see, 
e.g., Chiang & Choi 2008). 

4.3. Parent Body and Planet Origins 

Though the direct detection of parent bodies is beyond 
the reach of current observations, our study has provided 
some evidence that they reside mostly on nearly purely 
secularly forced orbits with small free eccentricities. In 
principle, they do not have to; we found by experimenta- 
tion in §3.3.3 that large free eccentricities are also com- 
patible with long-term stability if the particles are pro- 
tected by mean-motion resonances. 

How did the parent bodies choose one class of stable 
orbit over the other? The answer, as suggested also by 
Quillen & Faber (2006), likely involves collisional dissi- 


pation. Collisions dissipate random orbital motions and 
compel planetesimals to conform towards closed, non- 
intersecting orbits viewed in the frame rotating with the 
perturbation potential (e.g., Paczynski 1977; Goldreich & 
Tremaine 1982, section 5.4). If that perturbation poten- 
tial arises from Fom b, the special closed orbits include 
the secularly forced orbits that we have been highlight- 
ing throughout our study, but they do not include the 
resonant orbits that emerged from our experiment. Col- 
lisional dissipation and relaxation onto closed orbits oc- 
cur across the entire collisional cascade, up to the largest 
parent bodies (which by definition collide once over the 
system age). 

We do not expect the destructive nature of collisions to 
qualitatively alter this picture, since what is important 
here is the dissipation of random kinetic energy, and that 
occurs whether or not collisions are destructive. Post- 
collision fragments will have free eccentricities that are 
small compared to forced eccentricities, insofar as post- 
collision fragment velocities (measured relative to the 
center of mass) are small compared to ef or ced^i?7? ~ 400 
m/s. At least in the collisional genesis of the Eunomia 
and Koronis asteroid families in our solar system’s main 
belt, ejection velocities of the largest post-collision rem- 
nants range from 4 to 90 m/s (Michel et al. 2001). 

It is also possible that relaxation occurred while the 
parent bodies were forming, when collisions were gentler 
and agglomerative. The process of relaxation onto forced 
orbits can be explored using fast numerical simulation 
techniques for inelastically colliding, indestructible par- 
ticles (Lithwick & Chiang 2007); in fact, we have started 
to run such simulations and clearly observe relaxation. 

What are the origins of Fom b and the belt? It seems 
likely that the belt is what remains of the original disk 
material that went into building Fom b. If we take the 
minimum parent body mass of (estimated in §2.6) 
and augment it by a factor of 10 2 to bring it to cosmic 
composition, then the minimum primordial mass for the 
belt is ~1Mj. This is comparable to our upper mass 
limit for Fom b. A working hypothesis is that Fom b ac- 
creted in situ from a primordial disk of gas and dust; that 
the hydrogen gas of the original disk has either accreted 
into planets or photoevaporated; and that today the re- 
maining solids in the belt are grinding down to dust, the 
in-plane velocity dispersion of parent bodies excited so 
strongly by Fom b, and possibly other perturbers, that 
collisions are destructive rather than agglomerative. The 
last tenet is supported by our Figure 6, which shows free 
eccentricity dispersions that imply relative parent body 
velocities upwards of ~100 m/s (see also our §2.3). 
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